Effects of Neighborhood Environmental Circumstances on Thalamocortical Connectivity
Prepare data
Glasser parcel names
#Glasser regions with corresponding labels and tract names
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) S-A axis
S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.regions, by = "orig_parcelname")
S.A.axis$SA.axis.bin <- as.numeric(cut2(S.A.axis$SA.axis, g = 5)) #bins for quintile enrichment analysis
SA.colorvector <- c("#FFC228", "#FFD16A", "#e9dcf2", "#BE93C4", "#6F1382") #colors to use for plotting for each of the 5 quintiles
S.A.axis$SA.axis.quartile <- as.numeric(cut2(S.A.axis$SA.axis, g = 4)) #quartiles for sensitivity enrichment
S.A.axis$SA.axis.decile <- as.numeric(cut2(S.A.axis$SA.axis, g = 10)) #deciles for sensitivity enrichmentEnvironment statistics
setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/") #gam output dir, from gam_functions/fit_envGams.R
#list all .csv files in the environment results directory
files <- list.files(getwd(), pattern = '.csv', ignore.case = T, full.names = F)
for(i in 1:length(files)){
filename <- files[i]
Rname <- gsub('.{4}$', '', filename)
x <- read.csv(filename, header = TRUE)
assign(Rname, x)
}
rm(x)Spatial permutation (spin) test parcel rotation matrix
perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinningDataset-specific tract lists
#generated by tract_measures/tractlists/thalamocortical_tractlists.R
tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv")
tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv") Associations Between Neighborhood and Household SES and Thalamocortical Connectivity
Number of significant environmental effects
##Function to calculate the number and percent of thalamocortical connections showing a significant association with an environmental characteristic
sig.effects <- function(df, env.measure){
enveffects.df <- df
enveffects.df$significant <- p.adjust(enveffects.df$GAM.cov.pvalue, method = c("fdr")) < 0.05
sigeffect.totaln <- sum(enveffects.df$significant)
sigeffect.percent <- round(sigeffect.totaln/length(enveffects.df$significant)*100, 2)
print(sprintf("%s thalamocortical connections (%s percent) show a significant association between thalamocortical FA and the %s", sigeffect.totaln, sigeffect.percent, env.measure))
}PNC
Significant neighborhood-level environment effects (neighborhood SES)
sig.effects(df = envSES_maineffects_FA_glasser_pnc, env.measure = "neighborhood environment (envSES factor)")## [1] "128 thalamocortical connections (57.4 percent) show a significant association between thalamocortical FA and the neighborhood environment (envSES factor)"
controlling for parental education
sig.effects(df = envSES_maineffects_educontrolled_FA_glasser_pnc, env.measure = "neighborhood environment (envSES factor)")## [1] "125 thalamocortical connections (56.05 percent) show a significant association between thalamocortical FA and the neighborhood environment (envSES factor)"
Percent of positive significant neighborhood-level environment effects
(envSES_maineffects_FA_glasser_pnc %>% mutate(significant = p.adjust(GAM.cov.pvalue, method = c("fdr")) < 0.05) %>% filter(significant == TRUE) %>% filter (GAM.cov.tvalue > 0) %>% nrow() / 128) * 100## [1] 92.96875
Brain map of neighborhood-level environment effects
envSES_maineffects_FA_glasser_pnc <- envSES_maineffects_FA_glasser_pnc %>%
mutate(significant = p.adjust(envSES_maineffects_FA_glasser_pnc$GAM.cov.pvalue, method = c("fdr")) < 0.05)
sigeffects.df <- envSES_maineffects_FA_glasser_pnc %>% filter(significant == TRUE)
sigeffects.positivet.df <- sigeffects.df %>%
filter(GAM.cov.tvalue > 0) %>%
left_join(glasser.regions, ., by = "tract") %>%
select(label, tract, GAM.cov.tvalue)
sigeffects.positivet.df$tractlist <- sigeffects.positivet.df$tract %in% tractlist.pnc$tract
ggplot() +
geom_brain(data = sigeffects.positivet.df, atlas = glasser, mapping = aes(fill = GAM.cov.tvalue, colour=I("black"), size=I(.05))) +
theme_void() +
scale_fill_gradientn(colours = c(alpha("#323280", 0.15), "#672975"), limits = c(2, 5.25), oob = squish, na.value = "white") +
new_scale_fill() +
geom_brain(data = sigeffects.positivet.df, atlas = glasser, mapping = aes(fill = tractlist, colour=I("black"), size=I(.05))) +
scale_fill_manual(values = c("gray89", alpha("white", 0)), na.value = "white") +
theme(
legend.text = element_text(size = 5, family = "Arial", color = c("black")),
legend.key.size = unit(1, 'mm'),
legend.title = element_text(size = 5, family = "Arial", color = c("black")))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/envSES_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.35, height = 4)Significant household-level environment effects (average parental years of education)
sig.effects(df = householdSES_maineffects_FA_glasser_pnc, env.measure = "household environment (average parental years of education)")## [1] "16 thalamocortical connections (7.17 percent) show a significant association between thalamocortical FA and the household environment (average parental years of education)"
controlling for envSES
sig.effects(df = householdSES_maineffects_envSEScontrolled_FA_glasser_pnc, env.measure = "household environment (average parental years of education)")## [1] "0 thalamocortical connections (0 percent) show a significant association between thalamocortical FA and the household environment (average parental years of education)"
HCPD
Significant household-level environment effects (average parental years of education)
sig.effects(df = householdSES_maineffects_linear_FA_glasser_hcpd, env.measure = "household environment (average parental years of education)")## [1] "0 thalamocortical connections (0 percent) show a significant association between thalamocortical FA and the household environment (average parental years of education)"
Significant household-level environment effects (state-adjusted income-to-needs)
## [1] "0 thalamocortical connections (0 percent) show a significant association between thalamocortical FA and the household environment (INR)"
Thalamocortical Developmental Trajectories Stratified by Neighborhood Environment Factor Score
Trajectories in S-A axis quintiles
#Function to plot predicted developmental trajectories for 10th and 90th percentile envSES factor scores for bins of the S-A axis
plot_trajectories_byenv_SAaxis <- function(quintile, break.points, y.limits){
# Compute average developmental trajectories for low and high factors scores within each S-A axis quintile
df <- devtrajectories_byenvSES_FA_glasser_pnc %>% merge(., S.A.axis, by = c("tract"))
df$age <- round(df$age, 3)
df$envSES <- as.factor(df$envSES)
df <- envSES_maineffects_FA_glasser_pnc %>% select(tract, significant) %>% merge(., df, by = "tract")
df <- df %>% filter(significant == TRUE)
SAaxismooths.by.envSES <- df %>% group_by(envSES, SA.axis.bin, age) %>%
do(est.mean = mean(.$fitted)) %>%
unnest(cols = c(est.mean))
# Plot the trajectories for the requested quintile
SAtrajectories.byenv.plot <- SAaxismooths.by.envSES %>% filter(SA.axis.bin == quintile) %>%
ggplot(aes(x = age, y = est.mean, group = interaction(envSES, SA.axis.bin), linetype = envSES)) +
geom_line(color = SA.colorvector[quintile], linewidth = 0.4) +
theme_classic() +
labs(x = "\nAge (years)", y = "\n") +
scale_y_continuous(breaks = break.points, limits = y.limits) +
theme(
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))
return(SAtrajectories.byenv.plot)
}left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 1) %>%
mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
ggplot() +
geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"), size=I(.05))) +
theme_void() +
scale_fill_manual(values = SA.colorvector[1], na.value = "white")## merging atlas and data by 'label'
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile1.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 1, break.points = c(0.41, 0.42, 0.43, 0.44), y.limits = c(0.405, 0.44))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile1.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 2) %>%
mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
ggplot() +
geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"), size=I(.05))) +
theme_void() +
scale_fill_manual(values = SA.colorvector[2], na.value = "white")## merging atlas and data by 'label'
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile2.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 2, break.points = c(0.37, 0.38, 0.39, 0.40), y.limits = c(0.365, 0.40))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile2.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 3) %>%
mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
ggplot() +
geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"), size=I(.05))) +
theme_void() +
scale_fill_manual(values = "#d4bfe3", na.value = "white")## merging atlas and data by 'label'
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile3.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 3, break.points = c(0.33, 0.34, 0.35, 0.36), y.limits = c(0.325, 0.36))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile3.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 4) %>%
mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
ggplot() +
geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"), size=I(.05))) +
theme_void() +
scale_fill_manual(values = SA.colorvector[4], na.value = "white")## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry tract
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L_10pp EMPTY thalamus_L_10pp_autot…
## # ℹ 5 more variables: orig_parcelname <chr>, SA.axis <dbl>, SA.axis.bin <fct>,
## # SA.axis.quartile <dbl>, SA.axis.decile <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label tract orig_parcelname SA.axis SA.axis.bin
## 396 lh_L_10pp thalamus_L_10pp_autotrack L_10pp_ROI 282 4
## SA.axis.quartile SA.axis.decile
## 396 4 8
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile4.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry tract
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L_10pp EMPTY thalamus_L_10pp_autot…
## # ℹ 5 more variables: orig_parcelname <chr>, SA.axis <dbl>, SA.axis.bin <fct>,
## # SA.axis.quartile <dbl>, SA.axis.decile <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label tract orig_parcelname SA.axis SA.axis.bin
## 396 lh_L_10pp thalamus_L_10pp_autotrack L_10pp_ROI 282 4
## SA.axis.quartile SA.axis.decile
## 396 4 8
plot_trajectories_byenv_SAaxis(quintile = 4, break.points = c(0.33, 0.34, 0.35, 0.36), y.limits = c(0.33, 0.36))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/Devtrajectory_byenvSES_quintile4.pdf", device = "pdf", dpi = 500, width = 1.27 , height = 1)left_join(tractlist.pnc, S.A.axis, by = "tract") %>% filter(SA.axis.bin == 5) %>%
mutate(SA.axis.bin = as.factor(SA.axis.bin)) %>%
ggplot() +
geom_brain(atlas = glasser, mapping = aes(fill = SA.axis.bin, colour=I("black"), size=I(.05))) +
theme_void() +
scale_fill_manual(values = SA.colorvector[5], na.value = "white")## merging atlas and data by 'label'
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/SAquintile5.pdf", device = "pdf", dpi = 500, width = 2.3 , height = 3)## merging atlas and data by 'label'
plot_trajectories_byenv_SAaxis(quintile = 5, break.points = c(0.34, 0.35, 0.36, 0.37), y.limits = c(0.34, 0.375))Neighborhood Environment Effects Vary Across the Sensorimotor-Association Axis
S-A axis correlation
Spearman’s rho
spin.data <- left_join(spin.parcels, sigeffects.df, by = "tract") #merged data for spin tests, with all 360 parcels in the expected right hemisphere --> left hemisphere order (matching glasser.coords)
spin.data <- merge(spin.data, S.A.axis, by = c("label", "orig_parcelname", "tract"), sort =F)
cor.test(spin.data$GAM.cov.tvalue, spin.data$SA.axis, method = c("spearman"))##
## Spearman's rank correlation rho
##
## data: spin.data$GAM.cov.tvalue and spin.data$SA.axis
## S = 244280, p-value = 0.0005847
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3010667
Spatial permutation (spin) based p-value
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$GAM.cov.tvalue, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0.0056
Correlation plot
ggplot() +
geom_point(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue, color = SA.axis), size = 0.5) +
geom_smooth(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue), method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_y_continuous(limits = c(.5, 10)) +
scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
theme_classic() +
labs(x="\nS-A axis", y="Environment effect (t)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) Specificity test: cortical and thalamic axis correlations
#Glasser parcel x,y,z coordinates
load(file = "/cbica/projects/thalamocortical_development/Maps/parcellations/surface/hcp_mmp1.0.rda") #RDA file from the brainGraph package https://github.com/cwatson/brainGraph/blob/master/data/hcp_mmp1.0.rda with glasser region x, y, and z MNI coordinates. regions are in lh --> rh order
hcp_mmp1.0$x.mni[1:180] <- hcp_mmp1.0$x.mni[1:180]*-1 #convert R-> L hemisphere coords into medial lateral coords
hcp_mmp1.0$z.mni <- hcp_mmp1.0$z.mni*-1
hcp_mmp1.0$label <- NA
hcp_mmp1.0$label[1:180] <- glasser.regions$label[181:360]
hcp_mmp1.0$label[181:360] <- glasser.regions$label[1:180]
hcp_mmp1.0 <- hcp_mmp1.0 %>% select(-hemi)#C-Mt files for defining the group-level core-matrix gradient
CPt.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractaverage_CPt_primary.csv") %>% select("tract", "label", "orig_parcelname", "mean_CPt")
CPt.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractaverage_CPt_primary.csv") %>% select("tract", "label", "orig_parcelname", "mean_CPt")
CPt.gradient <- merge(CPt.glasser.pnc, CPt.glasser.hcpd, by = c("label", "tract", "orig_parcelname"), all.x = TRUE, all.y = TRUE, suffixes = c("_pnc", "_hcpd"))
CPt.gradient$CPt <- CPt.gradient %>% select(mean_CPt_pnc, mean_CPt_hcpd) %>% rowMeans(na.rm = T)
CPt.gradient <- CPt.gradient %>% select(CPt, label)
hcp_mmp1.0 <- left_join(hcp_mmp1.0, CPt.gradient, by = "label")#Function to compute the correlation between an environment map and an anatomical axis (x,y,z) as well as test the significance of the difference between the correlation with the anatomical axis versus the S-A axis
compute_axis_correlation <- function(axis){
df <- merge(spin.data, hcp_mmp1.0, by = "label", sort = F)
completeobs.n <- sum(!is.na(df$GAM.cov.tvalue))
#S-A axis - environment correlation (as computed above) for comparison of two overlapping correlations based on dependent groups
SA.axis.cor <- cor(df$SA.axis, df$GAM.cov.tvalue, use = "complete.obs", method = c("spearman"))
#Anatomical axis - environment correlation
anatomical.axis.cor <- cor(df[,axis], df$GAM.cov.tvalue, use = "complete.obs", method = c("spearman")) #correlation between anatomical axis and env metric
anatomical.axis.pvalue <- perm.sphere.p(x = df[,axis], y = df$GAM.cov.tvalue, perm.id = perm.id.full, corr.type = "spearman") #spin based p-value for the correlation
#S-A axis - anatomical axis correlation
SA.anatomical.cor <- cor(df$SA.axis, df[,axis], use = "complete.obs", method = c("spearman"))
#Test for significant difference between two correlations based on dependent groups (i.e., correlations with an overlapping input)
cocor.result <- cocor.pvalue <- cocor.dep.groups.overlap(
r.jk = SA.axis.cor,
r.jh = anatomical.axis.cor,
r.kh = SA.anatomical.cor,
n = completeobs.n, alternative = "two.sided", test = "hittner2003")
cocor.pvalue <- cocor.result@hittner2003[["p.value"]]
comparison.results <- data.frame(axis = axis, axis.cor = anatomical.axis.cor, axis.cor.pvalue = anatomical.axis.pvalue, SA.cor = SA.axis.cor, cocor.pvalue = cocor.pvalue)
return(comparison.results)
}Anterior-posterior axis
## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 y.mni 0.06846274 0.38685 0.3010667 0.0005890574
Medial-lateral axis
## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 x.mni 0.1841123 0.10335 0.3010667 0.3319559
Dorsal-ventral axis
## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 z.mni 0.2117744 0.07605 0.3010667 0.3909002
Core-matrix gradient
## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 CPt -0.1265908 0.15965 0.3010667 4.3867e-07
Correct axis correlation ps
axis.cor.ps <- c(0.0056, 0.38685, 0.10335, 0.07605, 0.15965) #S-A, A-P, M-L, D-V, C-M
p.adjust(axis.cor.ps, method = c("fdr"))## [1] 0.0280000 0.3868500 0.1722500 0.1722500 0.1995625
Correct correlation comparison ps
anatomical.ps <- c(0.0005890574, 0.3319559, 0.3909002, 4.3867e-07) #A-P, M-L, D-V, C-M
p.adjust(anatomical.ps, method = c("fdr"))## [1] 1.178115e-03 3.909002e-01 3.909002e-01 1.754680e-06
#Put above results into df for plotting
axis.results <- data.frame(Axis = factor(), corr = double())
axis.results <- axis.results %>% add_row(Axis = "S-A", corr = 0.3010667)
axis.results <- axis.results %>% add_row(Axis = "A-P", corr = 0.06846274)
axis.results <- axis.results %>% add_row(Axis = "D-V", corr = 0.2117744)
axis.results <- axis.results %>% add_row(Axis = "M-L", corr = 0.1841123)
axis.results <- axis.results %>% add_row(Axis = "C-M", corr = -0.1265908)
axis.results$Axis <- factor(axis.results$Axis, ordered = TRUE, levels = c("S-A", "D-V", "M-L", "A-P", "C-M"))
ggplot(axis.results, aes(x = Axis, y = corr, fill = Axis)) +
geom_bar(stat='identity') +
labs(x="") +
labs(y="\nEnvironment Effect Correlation") +
theme_classic()+
scale_fill_manual(values=c(alpha("#323280", 0.5), "#672975", "#672975", "#672975", "#672975")) +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))Effect magnitude enrichment testing (S-A axis quintiles)
Mean environment effect t-value by S-A axis quintile (empirical)
#Calculate empirical enrichment statistic for each S-A axis bin
envSES_effects_SAaxisbins <- spin.data %>%
group_by(SA.axis.bin) %>%
do(mean.tvalue.significant = mean(.$GAM.cov.tvalue, na.rm = TRUE)) %>%
unnest(cols = c("mean.tvalue.significant"))## # A tibble: 5 × 2
## SA.axis.bin mean.tvalue.significant
## <dbl> <dbl>
## 1 1 2.70
## 2 2 3.60
## 3 3 2.79
## 4 4 3.91
## 5 5 4.57
Null distribution and enrichment testing by S-A axis quintile
#Calculate null distribution of t-values by bin based on spherical spatial permutations (spins)
## inputs to spatially permute (i.e., spin) with the pre-computed spatial permutation matrix
x <- spin.data$SA.axis.bin #cortical map 1 to spatially permute: S-A axis bins
y <- spin.data$GAM.cov.tvalue #cortical map 2 to spatially permute: significant envSES-thalamocortical FA t-values
perm.id <- perm.id.full
## spin the empirical data
nroi = dim(perm.id)[1] #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = y.perm = array(NA,dim=c(nroi,nperm)) #initialize
for (r in 1:nperm) {
for (i in 1:nroi) {
x.perm[i,r] = x[perm.id[i,r]] #spinning x, spatially permuted bins
y.perm[i,r] = y[perm.id[i,r]] #spinning y, spatially permuted t-values
}
}
## compute mean t-value for each S-A axis bin using spatially permuted y inputs
### set up spun (y) and empirical (x) df
y.spatialperm <- as.data.frame(y.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(y.perm))))
y.spatialperm.empiricalx <- cbind(x, y.spatialperm) %>% as.data.frame()
colnames(y.spatialperm.empiricalx)[1] <- c("empirical.x")
### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedy <- y.spatialperm.empiricalx %>% group_by(empirical.x) %>% #group by empirical S-A axis bins
dplyr::summarize(across(everything(), \(x) mean(x, na.rm = TRUE))) %>% #calculate mean t-value for permuted data
select(-empirical.x) %>% t() %>% as.data.frame() %>%
set_names(c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue"))#Function to plot the null distribution of environment effect t-values for a given S-A axis bin and test whether the true t-value is significantly smaller or greater in magnitude than the null (permutation based p)
enveffects_SAbin_enrichment <- function(bin, enrichment_type){
nperm = 10000
SAbin.tvalue.empirical <- envSES_effects_SAaxisbins %>% filter(SA.axis.bin == bin) %>% select(mean.tvalue.significant)
#Histogram
SA.colorvector <- c("#FFC228", "#FFD16A", "#e9dcf2", "#BE93C4", "#6F1382")
null.tvalue.histogram <- envSES_effects_SAaxisbins_permutedy %>% select(sprintf("SA.bin%s.tvalue", bin)) %>%
ggplot(aes(x = .[,1])) + geom_histogram(fill = c(SA.colorvector)[bin]) +
theme_classic() +
geom_vline(xintercept = SAbin.tvalue.empirical$mean.tvalue.significant, linewidth = 0.25) +
xlab("Environment effect (t-value)") +
ylab("Number of nulls") +
scale_x_continuous(breaks = c(2, 3, 4, 5), limits = c(1.25, 5.25)) +
theme(
legend.position = "none",
axis.text.x = element_text(size = 5, family = "Arial", color = c("black")),
axis.text.y = element_blank(),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line.x = element_line(linewidth = .2),
axis.line.y = element_blank(),
axis.ticks = element_line(linewidth = .2))
#Permutation-based p
if(enrichment_type == "smaller"){
perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
}
if(enrichment_type == "larger"){
perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
}
permutation.results <- list(null.tvalue.histogram, perm.p)
return(permutation.results)
}S-A quintile 1
## [[1]]
##
## [[2]]
## [1] 0.0921
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile1_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)S-A quintile 2
## [[1]]
##
## [[2]]
## [1] 0.5832
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile2_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)S-A quintile 3
## [[1]]
##
## [[2]]
## [1] 0.0687
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile3_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)S-A quintile 4
## [[1]]
##
## [[2]]
## [1] 0.1859
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile4_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)S-A quintile 5
## [[1]]
##
## [[2]]
## [1] 0.025
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure7/quintile5_enveffect_spins.pdf", device = "pdf", dpi = 500, width = .95, height =.9)Enrichment sensitivity test: deciles
#Calculate empirical enrichment statistic for each S-A axis bin
envSES_effects_SAaxisbins <- spin.data %>%
group_by(SA.axis.decile) %>%
do(mean.tvalue.significant = mean(.$GAM.cov.tvalue, na.rm = TRUE)) %>%
unnest(cols = c("mean.tvalue.significant"))#Calculate null distribution of t-values by bin based on spherical spatial permutations (spins)
## inputs to spatially permute (i.e., spin) with the pre-computed spatial permutation matrix
x <- spin.data$SA.axis.decile #cortical map 1 to spatially permute: S-A axis bins
y <- spin.data$GAM.cov.tvalue #cortical map 2 to spatially permute: significant envSES-thalamocortical FA t-values
perm.id <- perm.id.full
## spin the empirical data
nroi = dim(perm.id)[1] #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = y.perm = array(NA,dim=c(nroi,nperm)) #initialize
for (r in 1:nperm) {
for (i in 1:nroi) {
x.perm[i,r] = x[perm.id[i,r]] #spinning x, spatially permuted bins
y.perm[i,r] = y[perm.id[i,r]] #spinning y, spatially permuted t-values
}
}
## compute mean t-value for each S-A axis bin using spatially permuted y inputs
### set up spun (y) and empirical (x) df
y.spatialperm <- as.data.frame(y.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(y.perm))))
y.spatialperm.empiricalx <- cbind(x, y.spatialperm) %>% as.data.frame()
colnames(y.spatialperm.empiricalx)[1] <- c("empirical.x")
### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedy <- y.spatialperm.empiricalx %>% group_by(empirical.x) %>% #group by empirical S-A axis bins
dplyr::summarize(across(everything(), \(x) mean(x, na.rm = TRUE))) %>% #calculate mean t-value for permuted data
select(-empirical.x) %>% t() %>% as.data.frame() %>%
set_names(c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue", "SA.bin6.tvalue", "SA.bin7.tvalue", "SA.bin8.tvale", "SA.bin9.tvalue", "SA.bin10.tvalue"))enveffects_SAbin_enrichment <- function(bin, enrichment_type){
nperm = 10000
SAbin.tvalue.empirical <- envSES_effects_SAaxisbins %>% filter(SA.axis.decile == bin) %>% select(mean.tvalue.significant)
#Permutation-based p
if(enrichment_type == "smaller"){
perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
}
if(enrichment_type == "larger"){
perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
}
return(perm.p)
}S-A decile 10
## [1] 0.0085
Enrichment sensitivity test: quartiles
#Calculate empirical enrichment statistic for each S-A axis bin
envSES_effects_SAaxisbins <- spin.data %>%
group_by(SA.axis.quartile) %>%
do(mean.tvalue.significant = mean(.$GAM.cov.tvalue, na.rm = TRUE)) %>%
unnest(cols = c("mean.tvalue.significant"))#Calculate null distribution of t-values by bin based on spherical spatial permutations (spins)
## inputs to spatially permute (i.e., spin) with the pre-computed spatial permutation matrix
x <- spin.data$SA.axis.quartile #cortical map 1 to spatially permute: S-A axis bins
y <- spin.data$GAM.cov.tvalue #cortical map 2 to spatially permute: significant envSES-thalamocortical FA t-values
perm.id <- perm.id.full
## spin the empirical data
nroi = dim(perm.id)[1] #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = y.perm = array(NA,dim=c(nroi,nperm)) #initialize
for (r in 1:nperm) {
for (i in 1:nroi) {
x.perm[i,r] = x[perm.id[i,r]] #spinning x, spatially permuted bins
y.perm[i,r] = y[perm.id[i,r]] #spinning y, spatially permuted t-values
}
}
## compute mean t-value for each S-A axis bin using spatially permuted y inputs
### set up spun (y) and empirical (x) df
y.spatialperm <- as.data.frame(y.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(y.perm))))
y.spatialperm.empiricalx <- cbind(x, y.spatialperm) %>% as.data.frame()
colnames(y.spatialperm.empiricalx)[1] <- c("empirical.x")
### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedy <- y.spatialperm.empiricalx %>% group_by(empirical.x) %>% #group by empirical S-A axis bins
dplyr::summarize(across(everything(), \(x) mean(x, na.rm = TRUE))) %>% #calculate mean t-value for permuted data
select(-empirical.x) %>% t() %>% as.data.frame() %>%
set_names(c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue"))enveffects_SAbin_enrichment <- function(bin, enrichment_type){
nperm = 10000
SAbin.tvalue.empirical <- envSES_effects_SAaxisbins %>% filter(SA.axis.quartile == bin) %>% select(mean.tvalue.significant)
#Permutation-based p
if(enrichment_type == "smaller"){
perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
}
if(enrichment_type == "larger"){
perm.p <- sum(envSES_effects_SAaxisbins_permutedy[,bin] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
}
return(perm.p)
}S-A quartile 4
## [1] 0.0323